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Abstract 

This  paper  describes  the  use  of  genetic  algorithms  (GAs)  for  the  optimal  design  of  phononic  bandgaps  in  periodic  elastic 
two-phase  media.  In  particular,  we  link  a  GA  with  a  computational  finite  element  method  for  solving  the  acoustic  wave 
equation,  and  find  optimal  designs  for  both  metal-matrix  composite  systems  consisting  of  Ti/SiC,  and  H20-filled  porous 
ceramic  media,  by  maximizing  the  relative  acoustic  bandgap  for  these  media.  The  term  acoustic  here  implies  that,  for  sim¬ 
plicity,  only  dilatational  wave  propagation  is  considered,  although  this  is  not  an  essential  limitation  of  the  method.  The 
inclusion  material  is  found  to  have  a  lower  longitudinal  modulus  (and  lower  wave  speed)  than  the  surrounding  matrix 
material,  a  result  consistent  with  observations  that  stronger  scattering  is  observed  if  the  inclusion  material  has  a  lower  wave 
velocity  than  the  matrix  material. 

©  2005  Elsevier  Ltd.  All  rights  reserved. 
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1.  Introduction 

Over  the  past  decade,  a  great  deal  of  attention  has  been  focused  on  the  design  of  phononic  bandgap  mate¬ 
rials  for  possible  use  in  a  wide  array  of  technological  applications,  such  as  improved  transducers,  ultrasound 
and  polarization  filters,  acoustic  barriers  for  noise  reduction,  and  vibrationless  environments  for  high-preci¬ 
sion  mechanical  devices.  Phonons  represent  quanta  of  elastic  vibrational  energy  which  govern  many  of  the 
physical  properties  of  solid  continua  such  as  thermal  conductivity  and  sound  propagation.  In  periodic  elastic 
media,  phononic  bandgaps  represent  frequency  regions  where  propagating  elastic  waves  do  not  exist.  Both 
theoretical  and  experimental  investigations  into  the  existence  of  absolute  acoustic  bandgaps  abound  in  the  lit¬ 
erature.  An  excellent  review  that  outlines  the  history  of  the  development  of  the  so-called  “band  theory”  for 
electrons,  photons,  and  phonons  can  be  found  in  Kushwaha  (1996). 
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A  variety  of  computational  methods  have  been  successfully  used  for  the  theoretical  prediction  of  phononic 
bandgaps  in  one-,  two-,  and  three-dimensional  elastic  media  consisting  of  periodic  inclusions  imbedded  in  a 
homogeneous  matrix.  Some  of  these  include  the  plane-wave  expansion  (PWE)  method  (Vasseur  et  al., 
2001;  Hou  et  al.,  2003;  Li  et  al.,  2003;  Zhang  et  al.,  2003a, b;  Kee  et  al.,  2000;  Suzuki  and  Yu,  1998;  Goffaux 
and  Vigneron,  2001),  the  finite  difference  time  domain  (FDTD)  method  (Sigalas  and  Garcia,  2000;  Tanaka 
et  al.,  2000;  Tanaka  and  Tamura,  2002),  the  multiple  scattering  (MST)  method  (Psarobas  et  al.,  2000;  Kafesaki 
and  Economou,  1999),  and  traditional  finite  element  (FEM)  methods  (Mead,  1996;  Zhang  et  al.,  2003c; 
Axmann  and  Kuchment,  1999).  The  computational  methods  solve  the  wave  equation  by  first  defining  the  ini¬ 
tial  geometry  and  mechanical  properties  of  the  inclusions  and  matrix  material.  Through  parametric  variation 
of  the  inclusion  and  matrix  properties  and  geometry,  various  solutions  to  the  forward  problem  described 
above  can  be  compared  with  experimental  results  (Montero  de  Espinosa  et  al.,  1998;  Henderson  et  al., 
2001;  Sheng  et  al.,  2003;  Miyashita  et  al.,  2003;  Meseguer  et  al.,  1999;  Diez  et  al.,  2000). 

Despite  the  rapidly  growing  literature  on  bandgap  materials,  only  a  few  applications  of  these  materials 
appear  in  the  literature.  Indeed,  very  few  papers  use  a  systematic  design  approach  to  create  PnBG  structures 
at  all.  Most  of  the  work  currently  being  undertaken  either  involves  the  computational  analysis  of  PnBG  struc¬ 
tures,  or  uses  an  ad  hoc  cut-and-try  approach  to  the  design  of  such  structures.  Notable  exceptions  can  be 
found  in  the  book  by  Bendsoe  and  Sigmund  (2003)  in  a  chapter  devoted  to  wave  propagation  problems 
and  applications  to  the  automotive  industry  (noise  reduction),  crashworthiness,  and  biomechanics.  The  opti¬ 
mization  of  structures  subjected  to  time-harmonic  (i.e.,  sinusoidal)  loading  has  been  applied  to  the  systematic 
design  of  phononic  bandgaps  (Sigmund  and  Jensen,  2002,  2003;  Jensen  and  Sigmund,  2003).  In  their  work, 
finite  element  methods  are  used  in  conjunction  with  the  method  of  moving  asymptotes  (Svanberg,  1987)  to 
optimize  the  relative  bandgap  size  in  both  finite  and  infinite  two-dimensional  periodic  elastic  media. 

Design  problems  for  optimal  cavity  shape  in  infinite  media  (Cherkaev  et  al.,  1998),  two-phase  periodic  elas¬ 
tic  media  (Grabovsky  and  Kohn,  1995;  Vigdergauz,  1997),  and  multifunctional  composites  (Torquato  et  al., 
2003)  have  been  studied  extensively  for  media  under  static  boundary  conditions;  much  less  work  has  been  pub¬ 
lished  for  the  optimal  design  of  such  media  subjected  to  time-harmonic,  or  transient  loading.  This  paper 
describes  a  new  approach  to  the  design  of  PnBG  structures  using  genetic  algorithms  (Holland,  1975,  1992; 
Goldberg,  1989;  Weile  and  Michielssen,  1997).  We  anticipate  that  this  research  will  lead  to  the  systematic 
design  of  PnBG  materials  for  a  wide  variety  of  applications,  such  as  the  shock  isolation  of  electronic  compo¬ 
nents  in  vehicles  or  aircraft  subjected  to  dynamic  loading,  piezocomposite  sensors  or  actuators,  and  acoustic- 
shields  or  housings  for  engine  noise  reduction.  Like  the  work  reported  in  Jensen  and  Sigmund  (2003),  and 
Sigmund  and  Jensen  (2002,  2003),  the  analysis  will  (at  least  initially)  be  based  on  the  finite  element  method, 
but  unlike  that  work,  we  use  a  local  simplex  algorithm  to  determine  the  relative  bandgap  size  from  a  popu¬ 
lation  of  GA  optimized  PnBG  structures.  Indeed,  some  work  in  this  vein  has  already  been  accomplished, 
and  has  been  reported  in  Gazonas  et  al.  (2004). 

The  remainder  of  this  paper  will  describe  a  novel  computational  method  for  PnBG  design.  Specifically,  Sec¬ 
tion  2  describes  a  simple  method  for  PnBG  analysis  based  on  the  finite  element  method.  This  section  will  serve 
to  highlight  some  of  the  challenges  in  creating  an  efficient  PnBG  analyzer  for  use  in  an  optimization  scheme. 
Section  3  describes  how  PnBG  structures  may  be  synthesized  using  an  analysis  method  such  as  that  described 
in  Section  2.  Finally,  Section  4  describes  actual  periodic  elastic  PnBG  structures  that  have  been  computation¬ 
ally  designed. 


2.  Analysis  of  PnBG  structures 

The  analysis  of  PnBG  structures  begins  with  a  description  of  the  problem.  For  simplicity,  this  discussion  is 
limited  to  scalar  dilatational  waves  in  two  dimensions,  but  it  can  easily  be  extended  to  vector  waves  or  three 
dimensions.  Fig.  1  shows  the  unit  cell  of  a  two-dimensional  PnBG  structure.  The  structure  shown  is  composed 
of  two  homogeneous  materials,  but  this  may  also  be  generalized  to  include  more  materials  or  even  inhomo¬ 
geneous  materials,  and  is  not  a  limitation  of  the  analysis  presented  here.  The  periodic  cell  itself  is  denoted  by 
£2,  and  is  defined  by  the  edge  lengths  tx  and  t2,  and  the  acute  angle  a.  The  edge  of  the  unit  cell  will  be  denoted 
by  r,  and  it  is  composed  of  four  parts  as  shown  in  the  figure:  TR,  TL,  Tt,  and  TB.  More  specifically,  TR 
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- -  h  - - 

Fig.  1.  Unit  cell  of  a  two-dimensional  periodic  structure. 


denotes  the  right  side  of  the  boundary,  FL  denotes  the  left  side  of  the  boundary,  FT  is  the  top  of  the  boundary, 
and  rB  is  the  bottom  of  the  boundary. 

The  scalar  wave  equation  governing  dilatational  wave  propagation  in  homogeneous  elastic  media  in  the 
absence  of  body  forces  can  be  written  as  (Graff,  1991), 

(c2y2  =  0  forrefl,  (1) 

where  A(r)  is  the  volumetric  strain  at  r,  and  c  =  yJ(A  +  2p)/p  is  the  speed  of  a  dilatational  wave,  A  and  p  are 
the  Lame  parameters,  and  p  is  the  material  density.  The  time-harmonic  scalar  wave  equation  can  also  be  writ¬ 
ten  in  terms  of  the  pressure  using  the  relation  p( r)  =  KA(r)  =  —  |  tr<7,  where  K  =  A  +  is  the  bulk  modulus,  a 
is  the  Cauchy  stress,  and  co  is  the  angular  frequency  of  the  disturbance, 

^V2  +  -^Jp( r)  =0  for  r  e  Q.  (2) 


The  boundary  conditions,  however,  are  less  standard.  Specifically,  since  the  periodic  structure  is  assumed  to  be 
composed  of  lossless  materials,  the  disturbance  must  not  grow  or  shrink  as  it  makes  its  way  across  the  lattice. 
Thus,  defining  the  lattice  vectors  fi  =  t\k  and  fi  =  i2  cos  ax  +  l2  sin  ay,  and  the  propagation  constant  k  of  the 
wave,  the  quasi-periodic  boundary  conditions  may  be  expressed  as  follows: 


p{ r  +  fi)  =  pir)^  for  r  e  FL, 
p( r  +  fi)  =  pif)^2  for  r  e  FB, 


dp(r  +  fi) 0p(r 


dn 


dn 


dp(r  + 12)  


dn 


dn 


for  r  e  rL, 
for  r  e  rB. 


(3) 


These  are  known  as  the  Floquet  boundary  conditions  (Brillouin,  1946;  Lee  and  Yang,  1973;  Auld,  1990). 

Given  k,  Eqs.  (2)  and  (3)  define  an  eigenvalue  problem  for  the  allowed  frequencies  of  propagation  co  for 
that  k.  Solving  this  problem  for  all  k  results  in  the  set  of  all  allowed  frequencies  of  propagation.  The  bandgap 
is  the  complement  of  this  set  (Collin,  1991;  Joannopoulos  et  al.,  1995;  Kushwaha,  1996). 

Clearly,  scanning  an  infinitely  broad  range  of  k  values  is  an  insurmountable  task.  Fortunately,  the  period¬ 
icity  of  the  problem  also  limits  the  meaningful  values  of  k  (Collin,  1991;  Joannopoulos  et  al.,  1995;  Kushwaha, 
1996).  In  particular,  the  reciprocal  lattice  vectors  kx  and  k2  are  defined  by  the  relation 

k m  ‘  fi  271  Smm  (4) 


where  Smn  is  the  Kronecker  delta:  Smn  =  1  if  m  =  n  and  Smn  =  0  if  m  ^  n.  It  can  be  shown  that  if  co  is  an  eigen¬ 
value  of  the  equation  for  a  given  wavevector  k,  then  it  is  also  an  eigenvalue  of  the  equation  for  -k,  and  also 
for  any  wavevector  of  the  form  k  +  mki  +  nk2  for  any  integers  m  and  n. 
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In  the  rest  of  Section  2,  we  describe  how  this  problem  can  be  solved  by  the  finite  element  method,  and  how 
the  solutions  can  be  used  to  find  the  bandgap  of  a  given  PnBG  structure. 

2.1.  Finite  element  formulation  for  PnBG  analysis 

In  this  subsection,  we  describe  a  method  for  finding  the  allowable  frequencies  of  propagation  for  a  given 
wavevector  using  the  finite  element  method  (FEM)  (Zienkiewicz,  1977;  Hughes,  2000).  First,  Eqs.  (2)  and 
(3)  are  recast  as  a  variational  problem.  The  solution  to  the  continuum  problem  may  be  written  as  a  variational 
problem  in  space  coordinates  (see,  e.g.,  Kohn  et  al.,  1972), 

311  =  3  J  J  \S7p\2  -  a>2^j-  dr  =  0.  (5) 

To  ensure  that  only  functions  which  are  properly  normalized  and  satisfy  the  Floquet  boundary  condition  are 
admitted,  admissible  functions  p  are  subject  to  ma xp  =  1  and 

P{ r  +  li)  =  p{r)e’k  ''  for  r  e  TL, 
p(x  + 12)  =  p( r)e/k'2  for  r  e  rB. 

Next,  Q  is  meshed,  and  p  is  approximated  by  finite  elements.  While  any  type  of  finite  element  discretization  can 
be  applied  in  principle,  Fig.  2  shows  linear  elements  on  triangles  and  will  form  the  basis  for  this  discussion. 
(Thus,  neither  the  simple  mesh  nor  the  assumption  of  linear  elements  is  a  restriction  of  the  technique,  only 
a  basis  for  discussion.)  In  this  discretization,  the  unknowns  are  located  at  the  corners  of  the  triangles  and 
the  unknown  function  varies  linearly  within  the  element.  Unknowns  are  indicated  in  Fig.  2  by  circles,  with 
black  circles  indicating  internal  nodes,  white  circles  indicating  independent  boundary  nodes,  and  gray  circles 
indicating  nodes  with  values  dependent  on  the  independent  boundary  nodes  through  the  Floquet  boundary 
condition.  Assuming  that  there  are  N  internal  nodes,  M  independent  boundary  nodes,  and  P  dependent 
boundary  nodes,  p  can  be  approximated  by  finite  elements: 

N  M  P 

p{ r)  ~  Y  VibA)  +  Y  Wibi+N(  r)  +  yy^jzibi+NjrM{  r).  (7) 

i- 1  i=  1  i=  1 

In  Eq.  (7),  the  vh  and  zt  are  complex  unknown  coefficients,  and  the  bf r)  are  basis  functions.  Substituting 
this  approximation  into  Eq.  (5)  yields  a  new  minimization  problem  over  the  space  of  the  finite  element 
discretization: 


Fig.  2.  Meshed  unit  cell  showing  internal  nodes  in  black,  independent  boundary  nodes  in  white,  and  dependent  boundary  nodes  in  gray. 


G.A.  Gazonas  et  al.  /  International  Journal  of  Solids  and  Structures  43  (2006)  5851-5866 


5855 


r  A  A  A  1 

^Mlll  ^liv  ^UW 

BUU  BUy  Buw 

1 

1 

V 

< 

A*  A  A 

^UV  ^vv  ^vw 

—  CD2 

DUV  Dvv  Dvw 

1 

W 

A*  A*  A 

_  ^UW  ^Y  \v  ^WW  _ 

R*  R*  R 

_  u\v  Dv\v  Dww  _ 

J 

1 

Z 

where  v,  w,  and  z  are  simply  vector  lists  of  the  appropriate  unknowns. 

Eq.  (8)  contains  the  discretization  of  the  partial  differential  equation  alone.  To  incorporate  the  boundary  con¬ 
dition,  the  values  of  the  P  dependent  boundary  nodes  in  the  vector  z  are  expressed  in  terms  of  the  values  of  the  M 
independent  boundary  nodes  in  the  vector  w  using  Eq.  (6).  This  procedure  results  in  an  equation  of  the  form 


z  =  Dw, 


(9) 


where  D  is  a  matrix  with  only  one  element  per  column.  Eq.  (9)  can  be  substituted  into  Eq.  (8),  and  the  min¬ 
imization  can  be  performed.  This  process  yields  a  generalized  eigenvalue  problem  to  be  solved  for  the  un¬ 
knowns  v  and  w,  and  the  frequency  cd\ 


Auu  Auv  +  AUWD 

a;  +  D*A*W  ( Aw  +  AVWD  +  D*  A;w  +  D*AWWD) 


V 

W 

=  CD 


Bn 


Buv  +  BUWD 


B;  +  D*B;w  (Bvv  +  BVWD  +  D*B;w  +  D*BWWD) 


V 

W 

(10) 


Once  this  problem  has  been  solved  for  a  given  k,  all  of  the  allowed  frequencies  of  propagation  and  modal  dis¬ 
tributions  are  known  for  that  wavevector. 

The  efficient  solution  of  the  generalized  eigenvalue  problem  of  Eq.  (10)  represents  one  of  the  most  challeng¬ 
ing  aspects  of  this  work;  currently,  our  solution  of  these  types  of  equations  uses  the  ARPACK  package 
(Lehoucq  et  al.,  1998)  for  sparse  eigenvalue  problems.  While  ARPACK  has  the  capability  to  solve  sparse, 
complex,  generalized  eigenvalue  problems,  it  operates  much  more  efficiently  for  standard  eigenvalue  problems; 
for  a  moderately  sized  problem  of  400  unknowns,  ARPACK  takes  about  12  s  to  compute  the  two  smallest 
eigenvalues  and  eigenvectors.  Our  implementation  uses  the  fact  that  the  B  matrix  of  Eq.  (10)  (i.e.,  that  matrix 
that  is  multiplied  by  cd2)  is  always  positive  definite  and  Hermitian  (Jin,  2002).  This  allows  it  to  be  decomposed 
using  the  Cholesky  algorithm  (Kincaid  and  Cheney,  1991).  This  decomposition  can  then  be  used  to  convert 
the  generalized  eigenvalue  problem  to  a  standard  eigenvalue  problem;  since  the  original  matrix  is  sparse,  a 
special  Cholesky  routine  was  written  to  do  this.  Once  this  conversion  has  been  accomplished,  ARPACK 
can  be  used  to  solve  the  standard  eigenvalue  problem.  This  reduces  the  computation  time  from  12  s  per  wave- 
vector  to  1  s  per  wavevector  (Gazonas  et  al.,  2004). 


2.2.  Determination  of  the  bandgap 


While  the  method  of  the  previous  section  is  a  necessary  component  of  a  bandgap  computation,  alone  it  is 
insufficient.  In  particular,  bandgap  determination  requires  the  location  of  those  areas  of  the  cd  axis  in  which  no 
mode  can  propagate.  To  accomplish  this,  it  is  necessary  to  determine  the  largest  and  smallest  values  of  cd  in 
each  band  (Fig.  4).  This  can  be  done  with  any  optimization  method  in  principle;  our  implementation  currently 
uses  the  simplex  method.  Once  each  band’s  support  is  known,  bandgaps  can  be  located  by  taking  the  comple¬ 
ment  of  the  union  of  the  support  of  each  band. 

Two  observations  are  crucial  here:  first,  since  an  optimization  method  is  needed  to  find  the  edges  of  each  band, 
that  optimization  method  must  be  both  efficient  and  robust.  Efficiency  can  generally  be  achieved  by  using  a  local 
optimization  method;  certainly,  genetic  algorithms  or  simulated  annealing  are  not  appropriate  here.  To  achieve 
robustness,  it  is  possible  to  run  the  local  optimizer  several  times  with  different  initial  starting  points  to  ensure  the 
true  band  edge  is  located.  Moreover,  if  the  PnBG  structure  has  any  symmetries,  the  band  edges  are  often  (but  not 
always)  located  at  band  edges.  These  points  can  be  tested  individually  to  further  guarantee  robustness. 

Second,  in  designing  the  PnBG  structure,  the  physical  quantity  of  interest  is  not  the  size  of  the  bandgap 
itself,  but  the  size  of  the  bandgap  relative  to  its  center  frequency.  This  is  because  the  structure  can  always 
be  scaled  to  change  its  band  of  operation.  Moreover,  this  implies  that  the  largest  bandgap  will  be  located 
between  two  of  the  first  few  bands.  This  limits  the  search  to  these  bands  alone. 
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Fig.  4.  Band  structure  for  the  PnBG  structure  of  Fig.  3.  The  planar  axes  graph  the  components  of  the  wavevector,  and  the  vertical  axis 
graphs  allowable  values  of  c o  in  rad/s.  The  bandgap  exists  in  the  interval  co  e  [1.4426, 1.9418],  corresponding  to  a  relative  bandgap  <yz~Q>1 


of  29.8%. 


Fig.  3  shows  a  meshed  PnBG  structure.  The  unit  cell  of  the  structure  is  a  unit  square,  and  it  is  meshed  with 
800  right  triangles.  Material  selection  is  shown  in  Fig.  3  where  white  corresponds  to  a  material  with  density 
1  kg/m3  and  longitudinal  modulus  of  1  Pa,  and  gray  corresponds  to  a  material  of  density  1  kg/m3,  but  longi¬ 
tudinal  modulus  0.0769  Pa.  (These  numbers  were  chosen  simply  for  illustration,  and  because  they  make  dis¬ 
cussion  easier.)  The  gray  square  in  the  center  of  the  mesh  occupies  exactly  one  quarter  of  the  unit  cell.  Using 
the  FEM  and  optimization  techniques  described  in  this  Section  and  its  predecessor,  the  size  of  the  relative 
bandgap  was  determined  to  be  29.8%  as  shown  in  the  band  diagram  of  Fig.  4.  Finally,  Figs.  5  and  6  show 
the  relative  distribution  of  acoustic  energy  at  the  bandgap  edges.  Specifically,  Fig.  5  shows  the  distribution 
of  energy  at  the  top  of  the  bottom  band,  and  Fig.  6  shows  the  distribution  of  energy  at  the  bottom  of  the 
top  band. 

3.  Synthesis  of  PnBG  structures 

The  primary  goal  of  this  work  is  the  creation  of  an  automated  method  for  the  design  of  PnBG  structures. 
While  analysis  methods  such  as  those  described  in  Section  2  are  necessary  in  such  an  endeavor,  they  are 
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Fig.  5.  Energy  distribution  at  co  =  1.4426  rad/s.  The  amplitude  of  the  mechanical  disturbance  in  this  plot  can  be  arbitrarily  scaled;  the 
color  axes  are  not  labelled,  but  red  corresponds  to  high  energy  density,  and  blue  to  low  energy  density. 
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Fig.  6.  Energy  distribution  at  co  =  1.9418  rad/s. 

insufficient.  A  design  technique  is  also  required.  Such  synthesis  requires  an  optimization  technique  such  as  the 
genetic  algorithm  (Holland,  1975,  1992;  Goldberg,  1989;  Weile  and  Michielssen,  1997). 

GAs  are  based  on  Charles  Darwin’s  Theory  of  Descent  with  Modification  by  Natural  Selection,  and  have  a 
number  of  advantages  over  more  traditional  optimization  techniques  in  practical  design  applications.  Among 
these  advantages  are: 

Resistance  to  local  optima.  As  stochastically  based,  global  optimization  algorithms,  GAs  tend  to  find  global 
optima  or  strong  local  optima  for  complicated  optimization  problems.  This  is  important  in  any  practical 
design  problem  to  ensure  the  quality  of  the  designs  created  by  the  optimizer. 

Lack  of  reliance  on  derivatives.  Genetic  algorithms  have  no  need  to  compute  the  derivatives  of  objective 
functions  to  optimize  them.  In  fact,  no  estimate  of  derivatives  is  used  directly  or  indirectly  by  the  GA.  Not 
only  does  this  avoid  any  computationally  expensive  derivative  calculation,  but  also  ensures  that  functions  with 
exceptionally  complex  contours  can  be  optimized  easily  and  efficiently. 

Ability  to  optimize  functions  of  continuous  and  discrete  parameters.  Very  often,  practical  design  problems 
entail  the  selection  of  both  parameters  that  are  allowed  to  vary  continuously,  and  other  parameters  that  can¬ 
not.  For  instance,  while  most  geometrical  parameters  in  the  design  of  a  PnBG  structure  can  take  any  value  in 
an  interval,  material  choice  may  need  to  be  limited  to  a  database  for  practical  reasons.  GAs  can  accomplish 
the  optimization  of  such  mixed  parameter  problems  with  no  need  to  approximate  continuous  parameters  by 
discrete  parameters. 
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Ability  to  be  implemented  in  parallel  GAs  are  a  member  of  a  class  of  so-called  “embarrassingly  parallel” 
algorithms  that  are  more  naturally  implemented  on  parallel  computers  than  on  serial  computers.  This  can 
be  essential  in  the  optimization  of  computationally  complex  problems.  Indeed,  most  GA  implementations 
experience  near  linear  speedups  when  ported  to  parallel  architectures. 

The  remainder  of  this  section  will  describe  GAs  as  they  apply  to  the  problem  of  designing  PnBG  structures. 
Section  3.1  will  discuss  GAs  in  general,  to  frame  the  discussion  of  how  these  algorithms  can  be  used  to  opti¬ 
mize  PnBG  structures.  Section  4  will  then  demonstrate  the  technique  with  an  actual  design  example. 

3.1.  Genetic  algorithms 

As  mentioned  in  the  introduction  to  this  section,  GAs  are  optimization  algorithms  based  on  Darwinian  the¬ 
ory.  GAs  differ  in  two  crucial  ways  from  more  standard  algorithms:  first,  they  operate  on  a  coded  version  of 
the  design  parameters,  instead  of  directly  on  these  parameters  themselves.  In  GA  terminology,  a  complete 
encoded  description  of  a  design  is  called  a  chromosome.  Second,  instead  of  optimizing  one  chromosome  at 
a  time,  GAs  optimize  an  entire  population  of  chromosomes  in  a  single  optimization  run  (Goldberg,  1989). 

Before  beginning  GA  optimization,  an  encoding  scheme  is  chosen  for  the  chromosomes.  Many  different 
representations  are  possible,  but  the  most  common  representation  is  binary.  Generally,  a  binary  chromosome 
is  composed  of  several  genes ,  each  of  which  is  composed  of  several  bits.  For  example,  in  this  work,  a  chromo¬ 
some  is  a  pair  of  square  arrays  of  bits,  signifying  the  material  choice  for  each  triangle  in  a  fixed  mesh  Fig.  7. 
Thus,  the  first  square  of  bits  in  Fig.  7  represents  the  state  of  one  of  the  triangles  in  each  of  the  25  parallelo¬ 
grams  created  by  a  rectilinear  mesh  of  the  cell,  and  the  second  square  represents  the  state  of  the  other  triangle. 
With  this  structure  we  have  a  possible  250  =  1015  possible  designs  in  the  GA  solution  space.  More  complex 
encodings  involving  more  material  choices,  etc.  are  also  possible  but  are  not  discussed  here. 

The  GA  begins  by  initializing  a  population  of  Np  chromosomes  at  random.  These  chromosomes  are  then 
evaluated  by  the  objective  function,  which  returns  a  measure  of  chromosome  goodness  in  maximization  prob¬ 
lems,  or  chromosome  badness  in  minimization  problems.  This  objective  function  represents  the  quantity  opti¬ 
mized  by  the  GA.  In  the  case  of  phononic  bandgap  design  presented  here,  the  objective  function  is  a  measure 
of  goodness:  the  relative  bandgap  discussed  above. 

After  initialization,  the  population  is  iteratively  subjected  to  the  three  genetic  operators’,  selection,  cross¬ 
over,  and  mutation.  Each  iteration  of  these  operators  in  sequence  (including  the  preceding  or  consequent  eval¬ 
uation  step)  is  termed  a  generation.  Selection  is  the  genetic  operator  responsible  for  implementing  the 
ubiquitous  Darwinian  dictum  of  survival  of  the  fittest.  The  application  of  selection  results  in  a  new  population 
of  Np  chromosomes  that  has,  on  average,  better  objective  function  values  than  the  original  population.  Again, 
this  can  be  accomplished  in  many  ways.  One  the  most  popular  of  these  is  binary  tournament  selection ,  which 
selects  two  candidate  chromosomes  at  random  for  each  spot  in  the  new  population,  and  places  the  better  one 
in  the  new  population.  On  average,  the  new  population  so  created  will  have  two  copies  of  the  best  chromo¬ 
some  in  the  population,  one  copy  of  the  median  chromosome,  and  no  copies  of  the  worst  chromosome. 

After  the  application  of  selection,  the  crossover  operator  is  applied.  Crossover  works  to  hybridize  the  chro¬ 
mosomes  that  survive  selection  by  combining  their  genetic  material.  In  a  standard  binary  coded  GA,  the  chro¬ 
mosomes  that  survive  selection  are  paired  randomly,  and  subjected  to  hybridization  with  a  probability  pc.  If  a 


Fig.  7.  Binary  chromosome  structure. 
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given  pair  is  to  be  hybridized,  a  random  point  between  bits  of  the  chromosome  is  chosen,  and  the  two  chro¬ 
mosomes  swap  genetic  material  after  that  point.  In  the  current  implementation,  if  a  given  pair  of  chromo¬ 
somes  is  to  be  crossed,  a  random  line  (horizontal  or  vertical)  is  drawn  in  the  bit  arrays,  and  the 
chromosomes  swap  the  material  either  to  the  right  of  or  below  the  line.  In  this  way,  the  hybridization  respects 
the  geometry  of  the  problem  at  hand.  Pairs  of  chromosomes  that  are  not  chosen  for  crossover  are  simply  left 
alone.  This  procedure  results  in  yet  another  population  of  Np  chromosomes,  containing  the  same  genetic  mate¬ 
rial  as  the  previous  generation,  but  jumbled.  Crossover  is  the  primary  search  modality  of  the  GA;  it  combines 
genetic  material  that  is  already  known  to  be  good  (since,  by  construction,  it  survived  selection). 

Finally,  the  mutation  operator  is  applied  to  randomly  alter  the  chromosomes  and  prevent  premature  con¬ 
vergence  to  a  suboptimal  result.  In  our  implementation,  each  bit  in  each  chromosome  in  the  population  is  sim¬ 
ply  negated  with  a  probability  pm.  (Generally,  pm  <  0.01  to  avoid  too  much  damage  to  the  chromosome.)  After 
mutation,  the  population  is  reevaluated,  and  a  new  generation  begins.  The  procedure  is  continued  until  a 
design  goal  is  met,  no  improvement  is  noticed  in  the  population,  or  a  fixed  number  of  generations  has  passed. 
Fig.  8  shows  the  design  coded  by  a  chromosome  at  a  given  slot  in  the  population  as  the  generations  progress. 
Initially,  the  design  is  completely  random  (Fig.  8(a)).  As  the  optimization  progresses,  suboptimal  genetic 


Fig.  8.  Evolution  of  an  optimal  design  showing  the  chromosome  structure  in  successive  generations,  from  the  initial  in  (a),  to  the  final 
in  (d). 
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material  is  eliminated  from  the  population,  which  becomes  more  homogenous  and  approaches  a  solution 
(Fig.  8(b)  and  (c)).  Finally,  an  optimal  (or  very  nearly  optimal)  design  is  found  and  the  entire  population  coa¬ 
lesces  on  the  solution  (Fig.  8(d)). 

Note  that  the  GA  operators  themselves  only  involve  very  trivial  computations.  The  most  expensive  part  of 
the  entire  GA  design  loop  is  the  evaluation  of  the  objective  function  itself.  Fortunately,  the  objective  function 
evaluation  can  be  done  in  parallel  very  easily:  Since  the  evaluation  of  each  chromosome  in  the  population  is 
independent  of  all  of  the  other  chromosome  evaluations,  the  entire  process  may  be  parallelized.  If  the  popu¬ 
lation  is  composed  of  1000  chromosomes,  1000  computers  can  be  used  to  evaluate  them,  resulting  in  factor  of 
1000  speedup  relative  to  a  serial  implementation  if  all  of  the  objective  function  evaluations  take  the  same 
amount  of  time.  This  perfect  linear  speedup  is  almost  never  realized  in  practice  since  not  all  chromosomes  will 
actually  take  the  same  amount  of  time  to  evaluate,  and  the  GA  operators  themselves  take  some  amount  of 
time.  Nonetheless,  it  is  not  uncommon  for  genetic  algorithms  running  in  parallel  on  N  <Np  processors  to 
achieve  a  speedup  of  0.95 N. 

4.  Design  of  PnBG  structures  using  GAs 

The  simple  genetic  algorithm  described  in  the  previous  section  can  be  combined  with  the  analysis  methods 
described  in  Section  2  to  create  an  automated  method  for  the  design  of  PnBG  structures.  This  section  dem¬ 
onstrates  how  to  do  this  for  a  fairly  simple  two-dimensional  case  consisting  of  a  matrix  material  and  an  inclu¬ 
sion  material. 

4.1.  PnBG  design  of  metal-matrix  composite  media 

A  GA  was  created  to  optimize  the  relative  acoustic  bandgap  for  a  Ti/SiC  metal  matrix  composite  system, 
proposed  for  use  in  next-generation  jet  aircraft  engines.  The  structure  was  assumed  to  be  arrayed  on  a  hex¬ 
agonal  grid  (i.e.,  in  Fig.  1,  a  =  60°)  composed  of  98  equilateral  triangles.  The  chromosome  was  composed 
of  two  7x7  arrays  of  bits,  with  each  bit  representing  the  material  choice  for  each  triangle.  With  this  structure 
we  have  a  possible  298  =  103°  possible  designs  in  the  GA  solution  space;  in  particular,  0  was  mapped  to  SiC, 
and  1  was  mapped  to  Ti.  The  GA  was  run  for  200  generations  using  a  population  size  Np  =  200,  and  required 
about  four  CPU-hours  on  64-nodes  of  a  parallel  Linux  cluster.  Note  that  the  results  presented  can  be  scaled; 
that  is,  the  unit  cell  can  shrink  by  a  factor  of  10  in  all  dimensions,  resulting  in  a  factor  of  10  increase  in  the 
frequencies  in  the  bandgap. 


Fig.  9.  A  genetic  algorithm  optimized  PnBG  structure  showing  the  final  disposition  of  the  Ti  inclusion  in  gray,  and  SiC  matrix  in  white. 
The  FE  analysis  mesh  discretization  is  16  times  finer  than  that  shown  in  the  figure. 
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The  structure  resulting  from  the  optimization  process  is  shown  in  Fig.  9.  It  has  a  relative  bandgap  of  8.15% 
as  shown  in  the  band  structure  diagram  of  Fig.  10. 

Finally,  the  distribution  of  energy  at  the  top  of  the  bottom  band  and  the  bottom  of  the  top  band  are  shown 
in  Figs.  11  and  12. 

To  ensure  that  these  results  were  truly  optimal,  the  GA  was  run  ten  more  times.  All  runs  returned  the  same 
optimal  relative  bandgap.  The  reason  for  this  robustness  is  made  clear  by  Fig.  13  which  shows  one  of  the  struc¬ 
tures  returned  by  these  runs.  A  careful  examination  shows  this  structure  to  be  identical  to  that  of  Fig.  9  aside 
from  a  shift  of  the  unit  cell,  and  a  reflection  over  the  x-axis.  In  short,  all  of  the  designs  returned  were  struc¬ 
turally  identical  to  the  original  design  with  the  material  with  the  lower  longitudinal  wave  speed,  Ti,  placed 
within  the  inclusion,  a  result  that  is  consistent  with  other  findings  (Economou  and  Sigalas,  1993;  Sigalas 
and  Economou,  1994;  Kafesaki  and  Economou,  1999).  It  is  interesting  to  note  that  the  design  which  maxi¬ 
mizes  the  relative  bandgap  for  this  material  system  consists  of  Ti  inclusions,  with  a  filling  fraction  of 
24.5%,  within  a  SiC  matrix,  and  is  antithetical  to  that  proposed  for  next-generation  jet  engine  design  which 
is  driven  by  strength  and  thermomechanical  performance  considerations  (Meguid  et  al.,  2002). 


x  104 


/fy  (rad/m)  kx  (rad/m) 

Fig.  10.  The  band  structure  of  the  GA  designed  PnBG  structure. 
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Fig.  11.  Energy  distribution  at  the  top  of  the  bottom  band. 
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Fig.  12.  Energy  distribution  at  the  bottom  of  the  top  band. 
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4.2.  PnBG  design  of  fluid-filled  porous  media 

Biot  (1956)  developed  a  theory  of  wave  propagation  in  fluid-saturated  porous  media  that  has  found  wide¬ 
spread  use  in  the  oil-industry  and  provides  the  framework  for  both  observational  (Plona,  1980),  and  more 
recent  theoretical  work  (Berryman  and  Wang,  2001;  Kim  et  al.,  2003;  Liu  et  al.,  2005;  Rajagopal  and  Tao, 
2005).  Despite  the  fact  that  porous  media  are  commonly  used  as  sound  absorbers  (Sgard  et  al.,  2005;  Rossetti 
et  al.,  2005),  we  are  not  aware  of  any  research  which  applies  formal  optimization  methods  to  the  systematic 
design  of  such  media.  In  this  section,  we  show  how  the  GA  can  be  used  to  find  optimal  designs  of  periodic, 
fluid-filled  porous  media  consisting  of  a  SiC  matrix,  with  H20  filling  the  pore  space.  Instead  of  fixing  the  unit 
cell  geometric  parameters,  a,  and  i2  as  was  done  in  the  previous  section,  we  can  generalize  our  approach  by 
having  the  GA  search  for  the  optimal  unit  cell  geometric  parameters  i\  and  a  illustrated  in  Fig.  1;  this  is  done 
with  the  realization  that  the  optimal  cavity  shape  tends  to  mimic  the  shape  of  the  unit  cell.  With  this  approach, 
the  GA  found  that  the  largest  relative  bandgap,  nearly  113%,  was  for  square-shaped  unit  cells  (and  square¬ 
shaped  pores),  t\  —  t2  =  1  m,  with  pore-filling  fractions  of  16%;  Li  et  al.  (2003)  also  find  large  bandgaps  for 
square,  rigid  rods  in  an  air  host  by  simply  rotating  the  axes  of  the  rods.  The  contour  plot  of  the  relative 


Fig.  15.  Optimized  fluid-filled  porous  medium  with  pore  in  gray  at  a  =  45°. 


Fig.  16.  Optimized  fluid-filled  porous  medium  with  pore  in  gray  at  a  =  90°. 
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bandgap  as  a  function  of  t\  versus  a  reveals  that  the  solution  for  the  shape  of  the  optimal  pore  occupies  an 
extremely  small  region  in  the  lower  right-hand  region  in  the  parameter  space  plot  (Fig.  14).  It  is  interesting 
to  note  that  the  GA  performed  well  in  finding  a  solution  in  this  parameter  space,  which  is  much  akin  to  finding 
a  needle  in  a  haystack,  and  is  a  daunting  task  for  any  global  optimization  method.  Examples  that  illustrate  the 
optimal  designs  for  pores  with  unit  cells  oriented  at  45°  and  90°  can  be  seen  in  Figs.  15  and  16,  respectively. 

5.  Summary 

This  paper  has  described  the  use  of  GAs  for  the  optimal  design  of  phononic  bandgaps  in  periodic  elastic 
media.  A  GA  was  combined  with  a  computational  finite  element  method  for  solving  the  acoustic  wave  equa¬ 
tion,  and  optimal  designs  were  created  for  both  metal-matrix  composite  systems  consisting  of  Ti/SiC,  and  for 
H20-filled  porous  ceramic  media  by  maximizing  the  relative  acoustic  bandgap.  The  inclusion  material  always 
had  a  lower  longitudinal  modulus  (and  lower  wave  speed)  than  the  surrounding  matrix  material,  a  result  con¬ 
sistent  with  observations  that  stronger  scattering  is  observed  if  the  inclusion  material  has  a  lower  wave  velocity 
than  the  matrix  material.  Interestingly,  the  shape  of  the  inclusions  were  always  regular  polygons  that  mim¬ 
icked  the  geometry  of  the  unit  cell  itself.  Because  this  result  is  likely  the  result  of  our  technique  for  describing 
the  inclusions  with  the  GA,  this  limited  the  generality  of  the  results.  In  the  future,  we  plan  to  generalize  the 
geometric  description  of  the  unit  cell  and  inclusion,  so  that  truly  optimal  designs  can  be  realized.  Generalizing 
the  computational  framework  to  include  multiphase  composite  materials  may  also  eliminate  the  requirement 
for  periodic  structures,  as  the  mechanisms  for  bandgap  formation  in  multiphase  materials  will  permit  the 
design  of  locally  resonant  materials  discovered  by  Liu  et  al.  (2000).  It  should  be  noted  that  even  in  this 
instance,  the  assumption  of  a  periodic  structure  does  not  limit  the  ability  of  the  code  to  design  locally  resonant 
structures.  Finally,  the  natural  extension  of  the  methodology  to  include  full  solid  mechanics  and  three  space 
dimensions  will  enable  design  of  more  general  structures  of  the  type  described  in  Economou  and  Sigalas 
(1993). 
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1  VIRGINIA  POLYTECHNIC  INST 
COLLEGE  OF  ENGRG 
RBATRA 

BLACKSBURG  VA  24061-0219 

1  COMPUTATIONAL  MECHANICS 
CONSULTANTS 
J  A  ZUKAS 
POBOX  11314 
BALTIMORE  MD  21239-0314 

1  KAMAN  SCIENCES  CORP 
D  L  JONES 

2560  HUNTINGTON  AVE  STE  200 
ALEXANDRIA  VA  22303 

1  APPLIED  RSRCH  ASSOC 
D  E  GRADY 

4300  SAN  MATEO  BLVD  NE 
STE  A220 

ALBUQUERQUE  NM  871 10 

1  INTL  RSRCH  ASSOC  INC 
D  L  ORPHAL 
4450  BLACK  AVE 
PLEASANTON  CA  94566 

1  AKT  MISSION  RSRCH  CORP 
M  EL  RAHEB 
23052  ALCALDE  DR 
LAGUNA  HILLS  CA  92653 

1  WASHINGTON  ST  UNIV 
SCHOOL  OF  MECH  &  MATL  ENGRG 
J  L  DING 

PULLMAN  WA  99164-2920 

2  WASHINGTON  ST  UNIV 
INST  OF  SHOCK  PHYSICS 
Y  M  GUPTA 

J  ASAY 

PULLMAN  WA  99164-2814 

1  ARIZONA  STATE  UNIV 
MECH  &  ARSPC  ENGRG 
D  KRAJCINOVIC 
TEMPE  AZ  85287-6106 


NO.  OF 

COPIES  ORGANIZATION 

1  UNIV  OF  DAYTON  RSRCH  INST 
NSBRAR 

300  COLLEGE  PARK 
MS  SPC  1911 
DAYTON  OH  45469 

1  DIRECTOR 
USARL 
K  WILSON 
FRENCH  DEA  1396 
ADELPHI  MD  20783-1 197 

1  TEXAS  A&M  UNIV 

DEPT  OF  GEOPHYSICS 
T  GANGI 

COLLEGE  STATION  TX  77843 

1  UNIV  OF  SAN  DIEGO 

DEPT  OF  MATH  &  CMPTR  SCI 
A  VELO 

5998  ALCALA  PARK 
SAN  DIEGO  CA  92110 

1  NATL  INST  OF  STAND  &  TECH 

BLDG  AND  FIRE  RSRCH  LAB 
J  MAIN 

100  BUREAU  DR  MS  861 1 
GAITHERSBURG  MD  20899-861 1 

1  BUCKNELL  UNIV 

DEPT  OF  MECH  ENGRG 
C  RANDOW 
LEWISBURG  PA  17837 

ABERDEEN  PROVING  GROUND 

73  DIR  USARL 

AMSRD  ARL  WM 
J  ALTHOUSE 
SKARNA 
J  MCCAULEY 
J  SMITH 
T  WRIGHT 
AMSRD  ARL  WM  B 
J  NEWILL 
M  ZOLTOSKI 
AMSRD  ARL  WM  BA 
AMSRD  ARL  WM  BC 
P  PLOSTINS 


NO.  OF 
COPIES 


ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


AMSRD  ARL  WM  BD 
P  CONROY 
B  FORCH 
RLIEB 

R  PESCE  RODRIGUEZ 
BRICE 

AMSRD  ARL  WM  BF 
D  WILKERSON 
AMSRD  ARL  WM  EG 
E  SCHMIDT 
AMSRD  ARL  WM  M 
S  MCKNIGHT 
AMSRD  ARL  WM  MA 
R  JENSEN 

M  VANLANDINGHAM 
E  WETZEL 
AMSRD  ARL  WM  MB 
M  BERMAN 
L  BURTON 
T  BOGETTI 
W  CHOWDHURY 
W  DE  ROSSET 
W  DRYSDALE 
A  FRYDMAN 
D  HOPKINS 
L  KECSKES 
T  H  LI 

M  MINNICINO 
B  POWERS 
J  TZENG 

AMSRD  ARL  WM  MC 
R  BOSSOLI 
S  CORNELISON 
M  MAHER 
W  SPURGEON 
AMSRD  ARL  WM  MD 
B  CHEESEMAN 
ECHIN 
B  DOOLEY 
C  FOUNTZOULAS 
G  GAZONAS 
J  LASALVIA 
P  PATEL 
J  SANDS 
CF  YEN 

AMSRD  ARL  WM  RP 
J  BORNSTEIN 
AMSRD  ARL  WM  SG 
T  ROSENBERGER 
AMSRD  ARL  WM  T 
AMSRD  ARL  WM  TA 
S  SCHOENFELD 


AMSRD  ARL  WM  TB 
P  BAKER 
R  SKAGGS 
J  STARKENBERG 
AMSRD  ARL  WM  TC 
R  COATES 
K  KIMSEY 
D  SCHEFFLER 
S  SCHRAML 
AMSRD  ARL  WM  TD 
S  BILYK 
T  BJERKE 
D  CASEM 
J  CLAYTON 
D  DANDEKAR 
M  GREENFIELD 
K  IYER 
BLOVE 

M  RAFTENBERG 
E  RAPACKI 
M  SCHEIDLER 
S  SEGLETES 
T  WEERASOORIYA 
AMSRD  ARL  WM  TE 
J  POWELL 
B  RINGERS 
G  THOMSON 


NO.  OF 

COPIES  ORGANIZATION 

1  DERA 

N  J  LYNCH 
WEAPONS  SYSTEMS 
BUILDING  A20 
DRA  FORT  HALSTEAD 
SEVENOAKS 
KENT  TN  147BP 
UNITED  KINGDOM 

1  ERNST  MACH  INTITUT 

H  NAHAME 
ECKERSTRASSE  4 
D  7800  FREIBURG  1  BR  791  4 
GERMANY 

1  FOA2 

P  LUNDBERG 
S  14725  TUMBA 
SWEDEN 

1  PCS  GROUP 

CAVENDISH  LABORATORY 
W  G  PROUD 
MADINGLEY  RD 
CAMBRIDGE 
UNITED  KINGDOM 

1  CENTRE  D  ETUDES  DE  GRAMAT 
J  Y  TRANCHET 
46500  GRAMAT 
FRANCE 

1  MINISTERE  DE  LA  DEFENSE 
DR  G  BRAULT 
DGA  DSP  STTC 
4  RUE  DE  LA  PORTE  DISSY 
75015  PARIS 
FRANCE 

1  SP ART  DIRECTION  BP  1 9 

DR  E  WARINGHAM 
10  PLACE  GEORGES 
CLEMENCEOUX 
9221 1  SAINT  CLOUD  CEDEX 
FRANCE 

1  LMT  CACHAN 
J  F  MOLINARI 

61  AVE  DU  PRESIDENT  WILSON 
94235  CACHAN  CEDEX 
FRANCE 


NO.  OF 

COPIES  ORGANIZATION 

1  TECHNICAL  UNIV  OF  CRETE 
G  EXADAKTYLOS 
DEPT  OF  MINERAL  RES  ENGNG 
CHANIA  CRETE 
GREECE 

1  ROYAL  MILITARY  COLLEGE  OF 
SCIENCE 

CRANEFIELD  UNIV 
J  MILLETT 

SHRIVENHAM  SWINDON 
SN6  8LA 

UNITED  KINGDOM 

1  UNIV  OF  MANCHESTER 
N  K  BOURNE 
PO  BOX  88 

SACKVILLE  STREET  MANCHESTER 
M60  1QD  UK 

1  BEN  GURIAN  UNIV  OF  NEGEV 
E  ZARETSKY 

DEPT  OF  MECH  ENG  BEER-SHEVA 
ISRAEL  84105 

2  RUSSIAN  ACADEMY  OF  SCIENCES 
INSTITUTE  FOR  HIGH  ENERGY 
DENSITIES 

G  I  KANEL 
S  V  RAZORENOV 
IVTAN  IZHORSKAYA  13/19 
MOSCOW  127412  RUSSIA 

1  INST  FUR 

M ATERI ALF ORSCHUN G  II 
C  TSAKMAKIS 
POSTFACH  3640 
FORSCHUNGSZENTRUM 
KARLSRUHE 
D  76021  KARLSRUHE 
GERMANY 

1  TECHNICAL  UNIV  OF  DENMARK 
DEPT  OF  MATHEMATICS 
M  BENDSOE 
2800  LYNGBY 
DENMARK 

1  TECHNICAL  UNIV  OF  DENMARK 
DEPT  OF  MECH  ENGRG 
O  SIGMUND 
2800  LYNGBY 
DENMARK 


NO.  OF 

COPIES  ORGANIZATION 


1  NATL  TECH  UNIV  OF  ATHENS 
DEPT  OF  ENG  SCI 
I  VARDOULAKIS 
ATHENS  15773 
GREECE 

1  UNIV  OF  PATRAS 

DEPT  OF  CIVIL  ENGRG 
D  BESKOS 
PATRAS  26500 
GREECE 

1  ARISTOTLE  UNIV 

OF  THESSALONIKI 
DEPT  OF  MECH  &  MATLS 
E  AIF ANTIS 
THESSALONIKI  54006 
GREECE 

1  DEMOCRITUS  UNIV  OF  THRACE 
DEPT  OF  CIVIL  ENGRG 
E  GDOUTOS 
XANTHI 
GREECE 


Intentionally  left  blank. 


